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Abstract 



Wc present an exact sampling method for the hrst passage event of a Levy process. 
The idea is to embed the process into another one whose hrst passage event can be 
sampled exactly, and then recover the part belonging to the former from the latter. The 
method is based on several distributional properties that appear to be new. We obtain 
general procedures to sample the first passage event of a subordinator across a regular 
non-increasing boundary, and that of a process with infinite Levy measure, bounded 
variation, and suitable drift across a constant level or interval. We give examples of 
, <~| application to a rather wide variety of Levy measures. 
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The first passage event of a Levyprocess is an important topic in applied probability and 
has received extensive study [cf. UH, l^-JJ, 21, 24, 36, and references therein]. However, 



■ although practically important and conceptually interesting, its exact sampling remains a 

£NJ . subtle issue. In particular, when the Levy measure of a process has infinite integral, except 

for a few well-known cases, the distribution of the process is not analytically available, 
which poses a significant hurdle to the exact sampling of its first passage event. 

Real-valued Levy processes with bounded variation form a large class. Since each such 
process is the difference between two independent subordinators, i.e., non-decreasing Levy 
processes, many properties of Levy processes with bounded variation can be deduced from 
those of subordinators. By themselves, subordinators not only play a significant role in 
the theory of Levy processes 0, U,[2^], but also have important applications 



As in the general case, for most subordinators with infinite Levy measure, exact sampling 
methods for the first passage event are lacking, although differential equations can be used 
to evaluate some parameters involved 
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In this paper, we present a method to sample the first passage event for a rather wide 
range of real-valued Levy processes with bounded variation. From now on, by sampling 
we always mean exact sampling. The method can sample 1) the first passage event of a 
subordinator across a regular non-increasing boundary, where the meaning of "regular" is 
specified in Section HI 2) the first passage event of a Levy process with bounded variation 
and non-positive drift across a positive constant level, and 3) the first exit event of a Levy 
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Figure 1: Sampling of the first passage event of a subordinator by embedding it into a 
"carrier" subordinator; see Section [1.11 for detail. 
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process with bounded variation and no drift out of a closed interval with as an internal 
point. As a by-product, the method also provides a way to sample infinitely divisible 
random variables alternative to one in 0]. 

We first give an overview of the method for subordinators, which is the main issue, and 
then an overview of its extension to real-valued Levy processes with bounded variation. 




1.1 Overview of method for subordinators 

In a nutshell, the idea is to embed a subordinator of interest into a "carrier" subordinator 
whose first passage event is tractable, and to utilize the latter to sample for the former. To 
put the idea into perspective, think of a Levy measure that can be decomposed as 

<p(x)dx + x{dx), with (p(x) = 1 {0 < x < r} ^e~ qx x~ a ~ x , (1.1) 

where r, 7 > 0, q > 0, a € (0, 1), and % ls a finite measure on (0, 00). This type of Levy 
measures coincide with those that can be decomposed as (p{x) dx + x{dx) with (p{x) = (7 + 
OixYjx -0-1 as x i 0, and give rise to many interesting processes discovered recently 




28|, |30|], such as Lamperti-stable process, whose Levy density is 1 {x > 0} e /3x (e x — 1) 
/3 < a + 1. For this particular process, we can set r = 00 in (jl.ip . However, in general, r 
has to be finite. 
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A subordinates with Levy density (jl.ip can be represented as the sum of two inde- 
pendent subordinators, one with Levy density ip, the other with Levy measure X- Since 
the latter is compound Poisson and only poses minor problems, we shall ignore it alto- 
gether here. Denoting by Z the subordinator with Levy density <p, it can be embedded 
into a stable subordinator as follows. Let X2 and X3 be independent subordinators which 
are also independent from Z and have Levy densities 1 {0 < x < r} (1 — e~ qx )x~ a ~ 1 and 
1 {x > r} x~ a , respectively. Then S = Z + X2 + A3 is a stable subordinator with index 



a 



35] . One can expect that the first passage event of S is tractable, which indeed is the 



case. Supposing the first passage event of S is sampled, we then have to recover the part 
due to Z from the sample values. This is the main issue we have to address. 

In general, let Z be a subordinator and c a non-increasing function on (0, 00). If Z has 
positive drift d > 0, then we can consider Z(t) — dt and c(t) — dt instead. Therefore, without 
loss of generality, let Z be a pure jump process. Our method requires some regularity of 
c. For now we ignore the issue and consider how to jointly sample the random variables 
T = {t > : Z(t) > c(t)}, Z(T—), and Z(T). Suppose the Levy measure of Z can be written 
as 1 {0 < x < r} e~ qx A(dx), where < r < 00 and q > 0, such that the Levy measure A(dx) 
gives rise to a subordinator S whose first passage event across any regular non-increasing 
boundary on (0, 00) can be sampled. Represent S as Z + X2 4- A3, such that Z, X2 and 
A3 are independent, with X2 and A3 having Levy measures 1 {0 < x < r} (1 — e~ qx )A(dx) 
and 1 {x > r}A(dx), respectively. Like Z, assume X2 and A3 are pure jump processes. 

The scheme of the method is shown in Fig. [TJ To start with, instead of c(t), let b(t) = 
c(t) A r be the "target boundary" for S to cross and r the corresponding first passage time. 
By assumption, we can sample (r, S(t— ), S(t)). Observe the following simple but crucial 
fact: since S(t—) < b(r) < r, S can only have jumps of size no greater than r in (0, r). 
Thus Z{t—) + X2(t— ) = S(t—). Now given t = t and S(t— ) = s < r, we have to sample 
Z(t—), which is possible for two reasons. First, the conditional distribution of Z(r— ) is the 
same as that of Z(t) given S(t) = s, as if t is fixed beforehand (cf. Section [3]). Second, using 
the properties of exponential tilting and upper truncation of Levy measures on (0, 00), the 
latter conditional distribution can be sampled (cf. Section [6]). In panel (a) of Fig. [TJ the 
sampled Z(r—) is less than S(t—). However, since X2 is compound Poisson, Z(t—) can be 
equal to S(r— ) with positive probability. Having got Z(r—), we next sample Z(r). The 
jump of S at r is A5 = S(t) — 5(r— ). By independence, only one of Z, X2 and X3 can 
have a jump at r, so Az is either or A$. Fig. [T] illustrates two scenarios. If As > r, then 
clearly it belongs to X3, giving Az = 0. If As < r, then by comparing the Levy measures 
of Z and X2, with probability e~ gAs (resp. 1 — e~ qAs ), As belongs to Z (resp. X2), giving 
Az = As (resp. Az = 0) (cf. Sections [3] and H|) . Thus Z{r) = Z(t—) + Az can be sampled. 
If Z(t) < c(r), then by strong Markov property, we can renew the procedure, but now with 
starting time point at r and starting value of S equal to Z(r). As can be expected, the 
procedure eventually stops, giving a sample value of (T, Z(T— ), Z(T)). 

Note that, if c is decreasing, then it is possible for S to creep across c, i.e., As = 0, 
as marked by * in panel (c). In this case, although Az = 0, if q > 0, it is possible that 
Z{t) < S(t) and so the procedure has to be renewed; see the scenario marked by B in 
panel (c). The characterization of creeping when c is linear is known 0, However, the 
case where c is non-linear appears to be still unresolved. 
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Figure 2: Sampling of the first passage event of Z = Z + — Z across a > 0, where Z^ are 
independent subordinators. 

a + Z~(r 3 ) 



a + Z~(r 2 ) ^ 



a + z-( n ) y 




z + (n) 



2 nd target boundary 

■ + (T2-) 




To implement the scheme, one can first sample r, then (S(t—),S(t)) conditional on 
r, and finally (Z(t—),Z(t)) conditional on (r, S(t— ), S(r)). Among these samplings, the 
one for r typically is the simplest. The other samplings require several theoretical results, 
which will be obtained in Section [3l Finally, it is easy to introduce a terminal point K < oo 
and sample (T", Z(T'—), Z(T')), with T 1 = T A K. In particular, if K = 1 and c = oo, the 
method samples an infinitely divisible random variable with upper truncated Levy measure 
1 {0 < x < r} e~ qx A(dx). As a passing remark, the so-called Vervaat perpetuity correspond 



to q = and A(dx) = adx, whose exact sampling is known id. Il5|. 



1.2 Overview of extensions 

The method described in Section 11.11 can be extended to Levy process with bounded vari- 
ation, as each such process is the difference of two independent subordinators. Fig. [2] illus- 
trates the sampling of the first passage event across a constant level a > by a Levy process 
with non-positive drift. Write the process as Z + — Z~ , where Z ± are independent subordi- 
nators, with Z + having no drift. The sampling can be thought of as having Z + to "catch 
up" with a + Z~ . To start with, let a be the target boundary for Z + to cross and t\ the cor- 
responding first passage time. It is evident that before n, Z^ stays below a + Z~. However, 
at ti, since Z + has a jump, it is possible for Z + to pass a + Z~ . We can use the method 
in Section [1.1 1 to sample Z + (t\). Meanwhile, since Z~ is independent of n, we can use the 
modification mentioned at the end of Section [TTT] to sample Z~(t\). If a+Z~ (n) < Z + (t±), 
then n is the first passage time of Z + — Z~ across a. If a + Z~(r\) > Z + (t\), then set 
a + Z~(ji) as the new target boundary for Z + to cross, this time with t\ as the starting 
time point and Z + (t{) the starting value for Z + . As long as lim^oo Z{t) = oo w.p. 1, the 
procedure eventually stops (cf. Section [5]). Here the assumption that Z has non-positive 
drift is critical, since otherwise Z + can creep across a + Z~ with positive probability. When 
this happens, the first passage times of Z + across the target boundaries shown in Fig. [2]will 
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Figure 3: A "phase plot" for the sampling of the first exit event of Z = Z + — Z out of a 
fixed interval [— a~,a + ], where > and Z^ are independent subordinators. 



converge but never be equal to the first passage time of Z across a, resulting in the iteration 
going on forever. It can be seen that if T = inf{t > : Z + (t) — Z~(t) > a}, then the proce- 
dure samples (T,Z + (T—),Z + (T),Z~(T)). As in Section fl.il a terminal point K < oo can 
be introduced so that one can sample (T", Z + (T'— ), Z + (T'), Z~(T')), where T' = T A K. 
If K < oo, the procedure eventually stops w.p. 1, whether or not lim^oo Z(t) = oo. 

It is also possible to sample the first exit event of Z = Z + — Z~ out of a fixed interval 
[— a~,a + ] with a ± > if it has no drift. To see how the method in Section 1 1.1 1 can be 
extended to this case, it is convenient to consider the "phase plot" of the Levy process 
W = (Z~ , Z + ), which shows the trajectory of W on the x-y plane without time axis. Then 
the first exit of Z out of [—a~, a + ] can be depicted as the first exit of W out of the band 
{(x,y) : —a~ < y — x < a + }. Suppose Z^ each satisfies the assumption in Section [l.ll so 
that, for example, Z + has Levy measure 1 {0 < x < r + } exp{— q + x}A + (dx), where A + (dx) 
is the Levy measure of a subordinator S + whose first passage event across any constant 
level can be sampled. To start with, let b~ = a" A r~ , b + = a + A r + , and set the top and 
right sides of the rectangle [0, b~] x [0, b + ] to be the target boundary. In panel (a) of Fig. [3l 
r~ > a~ and r + < a + , resulting in the rectangle as shown. Now sample the first passage 
event of S = (S~, S + ) across the target boundary. To do this, we can first independently 
sample the first passage times of across 6 ± . If, as shown in panel (a), S~ makes a 
crossing at time r before S + , then sample (S~(t— ), S~(t)), and sample S + (t—) = S + (t) 
conditional on «S + (r) < b + . We next can use the method in Section [1.1 1 to recover Z ± (r). 
In the scenario shown in Fig. [3l since the jump of S~ at r is greater than r~, it is not 
part of Z~ , and so we end up with W(r) as in panel (b). The procedure is then renewed. 
As long as Z ^ 0, the procedure will stop eventually and we get the first exit event of Z 
(cf. Sections [3] and [5]) . It can be seen that if T = inf{t > : Z(t) [— a~,a + ]}, then the 
procedure samples (T, Z ± (T— ), Z ± (T)). As in Section ll.l[ a terminal point K < oo can 
be introduced so that one can sample (T", Z ± (T / — ), Z ± (T')) with T' = T A K. 




(a) 



(b) 
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1.3 Organization of the paper 

Section [2] fixes notation and recalls preliminary results. Section [3] obtains distributional 
properties needed by the above methods. The basic tools are results on subordinators and 
fluctuation theory [cf. 0] , although some recent developments on the first passage event of 



a general Levy process could be exploited [cf. ]ll|, ]l2j, ]14], [29(. In addition, we need to get 
some detail on various conditional distributions of a subordinator as well as on creeping. 
Sections 0H6] present procedures to implement the methods and show their validity, and 
identify major sampling issues involved. Sections [7HH] show examples of application of the 
procedures to several types of Levy measures. The first type is given in (jl.ip . The second 
type consists of finite sums of Levy measures of the form (jl.ip . which require additional 
techniques. The third type consists of 1 {0 < x < r} x e~ x dx + x(dx), with x a finite 
measure on (0,oo). This type gives rise to the aforementioned Vervaat perpetuity and the 
Beta process in survival analysis, which has Levy density 1 {x > 0}e~ cx /(l — e~ x ) with 
c > and belongs to the Beta-class processes 3, 2^, 25, 31]. Rejection sampling and the 
Dirichlet distribution play an important role in these examples 

2 Preliminaries 

For x = (xi, . . . ,Xd) € M. d , denote by ||x|| its L 1 norm ^ \xi\. The convention inf = co 
will be used all along; "p.d.f." will stand for "probability density function with respect 
to Lebesgue measure". For a, b > 0, Beta(a, b) denotes the distribution with p.d.f. 
1{0 < x < ljx^^l - x) b ~ l /B(a,b), where B(a,b) = T(a)T(b)/T(a + 6), Gamma(a,6) 
denotes the one with p.d.f. 1 {x > Q}x a - l e- x / b /[b a T(a)\, and Exp(6) denotes the one with 
p.d.f. 1 {x > 0} e~ x l h jb. Finally, Unif(0, 1) denotes the uniform distribution on (0, 1). 

Let v and v be two probability measures on a measurable space f2 satisfying dv jdv oc 
g < C, where g > is a known function and C > a known constant. Then v can 
be sampled as follows: keep sampling £ ~ v and U ~ Unif(0, 1) until CU < This 
procedure is the well-known rejection sampling, which is exact and stops w.p. i B HQ. 



Let v be an infinitely divisible distribution on (0, oo) with Levy measure A. Given q > 0, 
Aq(dx) = e~ qx A(dx) is known as an exponentially tilted version of A. If v„ is the infinitely 
divisible distribution with Levy measure A q , then v q (dx) oc e~ qx v(dx) 0, 0, H, 0]- By 
setting g(x) = e~ qx l {x > 0}, rejection sampling can be used to sample v q based on v. 

The Dirichlet distribution Di(ai, . . . , a^) with parameters a\, . . . , a& > is a generaliza- 
tion of the Beta distribution. It can be defined as a distribution on M fc , such that for any 
measurable function g > on M fc and u ~ Di(oi, . . . , a^), 

EbMi = r r ( rt"rT Qfc ? / 1 ^ n xi ^ °> ^ n x T~ i ■ • • dx k~i, 

1 (ai) ■ ■ ■ 1 [a k ) J ^ 

where in the integral x^ = 1 — x% ^fc-i instead of a variate, and x = {x\ , . . . , x^) ■ The 

distribution has p.d.f. r(oi + • ■ -+a k ) nf=i x t'^ 1 1 11^=1 ^( a i) with respect to the degenerate 

measure a^idx) = 1 {xi > 0, ||x|| = 1} dxi ■ ■ ■ dx^_i 5{dxj, + x\ H h Xfc„i — 1), where 5 is 

the Dirac measure at 0. For convenience, we will refer to it as the p.d.f. of Di(oi, . . . ,a k )- 
Also, if k = 1, then for a > 0, define Di(a) to be the Dirac measure at 1. 
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3 Distributions of the first passage event 



We need several distributional properties to implement the method introduced in Section [TJ 
The main issue is, provided a subordinator Z can be embedded into another subordinator 
S = Z + X2 + -X3, how to recover the first passage event of Z from the one sampled for S. 
As noted in Section 11.11 we need to take into account the possibility that S creeps across a 
boundary. Also, for both the method for subordinators and its extension to Levy processes, 
we need to make sure the corresponding sampling procedures eventually stop. 

3.1 Results for subordinators 

Let X = (X\, . . . , Xd) be a Levy process taking values in [0,oo) d with Levy measure II 
and Laplace exponent f (1 — e~ ( ~ 9 ' x ^)H(dx), 9 G [0, oo) d . In the application later, X\ = Z, 
d = 3, and are independent. By definition, II has no mass at {0}. Let Ax be the jump 
process of X. The process S = \\X\\ = X\ + • • • + X^ is a subordinator with Levy measure 
IIs(d2;) = Jj v d 1 {||z|| G dx}II(dz), x G (0,oo), and jump process As = ||Ax||- In the 
rest of this section, we shall always assume 

n s (0,oo) = 00. (3.1) 

Denote Hs(x) = Hs( x , 00), and given c € C(0, 00), 

r c = rf = inf {t > : S{t) > c(t)} . (3.2) 

If c(t) = a > 0, the notation r a will be used instead of r c . To implement the method 
introduced in Section [LT| we will first sample r c , which is often easy, then (S(t c —), As(t c )) 
conditional on r c , and finally (X(t c — ), Ax{t c )) conditional on (r c , S(t c —), As(t c )). Among 
the results below, part 2) of Theorem 13.11 will be used for the conditional sampling of 
(X(t c — ), Ax(t c )), while Theorem 13.71 for the conditional sampling of (S(t c —),As(t c )). 

Theorem 3.1. Let c G C(0, 00) be non-increasing with c(0+) > and (I3.ip hold. 
1) LetQ = (0,oo) x [0,oo) d x ([0,oo) d \{0}). For (t,u,v) € 0, 

P{r c € dt, X(t c -) E du, A x {t c ) € dv} 

= 1{0 < c(t) - \\u\\ < \\v\\}dtP{X(t) G dn}n(dv), (3.3) 
P{r c G dt, X(r c -) G du, A x {r c ) = 0} 

= P {t c g dt, 5(t c ) = c(r c )} P {X(i) G du I 5(t) = c(t)} . (3.4) 

For t > 0, s G [0, c(t)], 2 > c(i) - s, and u, v G [0, oo) d , 

P {X(r c -) G du, A x (r c ) G dt» | r c = t, 5(r c -) = s, A s (t c ) = z) 

= P {X(t) G dn I S(t) = s} H z (dv), (3.5) 

with H z (dv) = P{V G dv | ||V|| = z}, where V is a random vector following the distribution 

v a {dv) = l{\\v\\ > a}n(du)/n s (a), 
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with a £ (0,-z) a fixed number. The conditional probability measure IL z (dv) is independent 
of the choice of a € (0, z). Furthermore, 



P {X(t c -) € du j r c = t, A s (r c ) = 0} = P {X(t) € du \ S{t) = c{t)} . (3.6) 

Proof. 1) It is clear that < r c < 00 w.p. 1. We first show f|3. 3|) . Following the proof for 
Proposition III. 2 in 0], let / > be a Borel function on such that f(t,u,v) = when 
\\v\\ = c(t) — \\u\\. Then 

/(r c ,I(r c -),A x (r c )) = £ f(t, X(t-), A x (t))l {0 < c(t) - S(t-) < \\A x (t)\\} . (3.7) 

t 

For each t > 0, define function H t (v) = f(t,X(t-),v)l{0 < c(t) - S(t-) < \\v\\} on 
[0, oo) d . Since H = (Ht) is a predictable process with respect to the filtration generated by 
Ax, by (|3.7p and the compensation formula, 

E[f(T c ,X(T c -),A x (r c ))] 

f'oo r r 

= J dtE J f(t,X(t-),v)l{0 < c(t) - S(t-) < \\v\\}U(dv) 

{ = } J dt J f(t,u,v)l{0 < c(t) - ||u|| < \\v\\}P{X(t) G du}U{dv) 

= / 1{0 < c(t) - ||u|| < \\v\\} f(t,u,v) dtP {X(t) € du}U(dv), 
Jo, 

where (a) is due to X(t—) ~ X(t) for any t > 0. Since / is arbitrary, this shows (j3.3|) for 
(t,u,v) € with II?; 1 1 ^ c(t) — \\u\\. It remains to consider (t,u,v) G O with ||«|| = c(t) — \\u\\. 
In this case, the right hand side of fj3.3[) is 0. If we define f(t, u,v) = 1 {v = c(t) — u > 0} 
for t > 0, n > and > 0, then by similar argument as above based on the compensation 
formula, but directly applied to S, 

P{S(t c -) < S(t c ) = c(r c )} = dt j P{S{t) € du}U s ({c(t) - u}). 

For each t, there is only a countable set of u with Us({c(t) — u}) > 0. On the other hand, 
under assumption (|3.ip . the distribution of S(t) is continuous, i.e., P{S(i) = u} = for any 



361 . Theorem 27.4]. As a result, J P{S(t) G du}Tls({c(t) — u}) = for each i, and so the 



multiple integral is 0. Finally, the proof of (|3.3p is complete by 

P {A x (t c ) + 0, 5(r c ) = c(r c )} = P {5(r c -) < 5(r c ) = c(r c )} = 0. (3.8) 

Now consider (|3.4p . Under f)3. 1 [) . 5 is strictly increasing w.p. 1. Clearly, A x (t c ) = 
implies S(t c ) = c(r c ). On the other hand, from (|3.8p . on the event iS(r c ) = c(r c ), Ax(t c ) = 
w.p. 1. Define r* = inf{t > : S(t) = c(t)}. Then w.p. 1, 

{r* < 00} = {t c = r*} = {S(r c ) = c(r c )}. (3.9) 
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Let / > be a Borel function on (0, oo) x [0, oo) d with bounded support. Then 
E[/(r c , X(t c — ))1 {S(t c ) = c(r c )}] can be expressed in two ways. First, from (|3.8p . it equals 

J f(t,u)l{\\u\\ = c(t)}P{T c £dt, I(r c -)edu, A x (r c ) = 0}. (3.10) 

Second, from (|3.8p and (|3.9p . the expectation also equals 

E[/(r c , X(r c ))l{,S(r c ) = c(r c )}] = E[/(r*, X(t*))1{t* < oo}] 

= y E[/(t, X(t)) | r* = t] P{n g dt} = J f(t, u) P{X(t) g du | r* = t}P{n g dt}. 

From the definition of r* and (|3.9p , the last integral is equal to 

J f(t, u) P {r c G dt, 5(t c ) = c(r c )} P G du | S{t) = c(t)} . (3.11) 

Since / is arbitrary, comparing the integrals in ()3.10|) and (|3.1ip then yields 

1 {||u|| = c(t)} P {r c G dt, X(t c -) G du, A x (t c ) = 0} 

= P {t c G dt, S(t c ) = c(r c )} P {X(t) G du | S(t) = c(t)} . 

Since the qualifier 1 {\\u\\ = c(t)} can be removed from the identity, (|3.4p follows. 

2) The process (X, S) is a Levy process with Hr x g\(dv, dz) = IL(dv) 5(dz — \\v\\). With 
similar argument as 1), for t > 0, s > 0, z > 0, u G [0, oo) d , and v G [0, oo) d \ {0}, 

P{r c G dt, X(r c -) G du, 5(t c -) G ds, A x (r c ) G du, A s (r c ) G dz} 

= 1 {0 < c(i) - s < z} dt P {X(i) G dn, S(t) G ds} n(du) 5(dz - 

On the other hand, applying (I3.3P directly to 5, we get 

P {t c G dt, S(r c -) G ds, A 5 (r c ) G dz} = 1 {0 < c(t) - s < z] dt P {5(t) G ds} n s (dz). 

Therefore, in order to get (j3.5[) . it suffices to show n(di>)<5(dz — ||u||) = n z (du )Hg(dz), which 
is equivalent to saying that for any measurable A C ([0, oo) d \ {0}) x (0, oo), 

j 1 {(v, z) G A} n(di>) 5{dz — \\v\\) = f 1 {(«,*) E A}U z (dv)U s {dz), 

The left hand side is J l{(f, \\v\\) G j4}n(df). To evaluate the right hand side, for 
any e > 0, let A e = {(v,z) G A : z > e}. Observe that 1 {z > e}Hs(dz)/Hs(s) is the 
distribution of ||V|| for V ~ v e . Then, from the definition of n z (du), 

J l{(v,z) eA £ }U z (dv)U s (dz) 

= n s (e) J 1 {(v, z) G A £ } P{V G dv | ||F|| = z}P{\\V\\ G dz} 

= S s (e) y 1 {(«, |H|) G P{y G du } = y 1 {(u, |H|) e A £ } 1 {\\v\\ > e} U(dv). 

Let e ! 0. Since j4 £ f A, the last integral converges to J" l{(u, ||f||) G A}n(di>). This 
completes the proof of (|3.5|) . Finally, since Aj(r c ) = if and only if As(t c ) = 0, and by 
(1331) . P{A s (r c ) / 0, 5(r c ) = c(r c )} = 0, (J3SD follows from <ET3]I . □ 



9 



Corollary 3.2. For a > and t > 0, define 

M*)= f n 5 (a-n)P{5(t) €du}. (3.12) 
Then, under the same assumption as Theorem \3.1\ 

r c(t) _ 

P{t c G dt, S{t c ) > c(r c )} = V c(t) (t) dt = dt / n s (c(t) - u) P {5(t) € du} . (3.13) 

J o 

/n particular, for constant a G (0,oo), T a /ias p.d/. ipa{t)- 
Proof. Apply <|3.3|> in Theorem 13.11 directly to 5 to get 

P{r c € dt, S(t c ) > c(r c )} = dt J 1{0< c(t) - u < v} P{S{t) G du}n s (d^), 

which is (I3J3D . Since P {S(t ) > a} = 1 @, Theorem III.4], P{r a G dt} = t/> a (t) dt. □ 



Definition 3.3. S is said to satisfy the continuous density condition, if S(t) has a p.d.f. gt 
on (0, oo) for each t > and the mapping (£, x) — > gt (x) is continuous on (0, oo) x (0, oo). 

We next obtain the p.d.f. of r c at the event that S creeps across a differentiable segment 



of c. For linear c, the result is shown in |17l ]. The following lemma is proved in Appendix. 

Lemma 3.4. Under the continuous density condition on S, the mapping (a,t) — > ip a {t) is 
continuous on (0, oo) x (0, oo), where tp a (t) is the p.d.f. of T a in (j3. 12[) . 

Proposition 3.5. Let c G C(0, oo) be non-increasing with c(0+) > 0. If c is differentiable 
on an open non-empty G C (0, oo) and S satisfies the continuous density condition, then 

P {r c G dt, 5(r c ) = c(r c )} = -c'(t)g t (c(t)) dt, t G G. (3.14) 

Proof. It suffices to consider t G G with c(t) > 0. Given such t, a := c(t) is fixed. Letting 
q(e) = P {t - e < t c < t}, for e > 0, we have q{e) = P {S(t - e) < c(t - e), 5(t) > a} = 
?i(e) + q 2 (s), where gi(e) = P {S(t - e) < a < S(t)}, q 2 (e) = P {a < 5(i - e) < c(t - e)}. 
Let r a = inf {s > : 5(s) > a}. Then ^(e) = P {t - e < r a < t}. By Corollary E21 and 
Lemma 13.41 ^ a is the continuous p.d.f. of r a , so the function P{r a < •} is differentiable 
with derivative ip a {t) at t. Then qi(e)/s —¥ ip a (t) = ifj c ^(t) as e 4 0. On the other hand, 

92(e) = Jq 6 ^ e ^ c( '^ <7(_ £ (a + x) dx. Since (t,x) — > gt(x) is continuous on (0, 00) x (0, 00) and 
c is differentiable at t, q 2 {e)je —¥ —c'{t)g t {c{t)) as e \. 0. We thus get 

toff = £ " = -d{t)gt(c(f)) +ip c(t )(t). 

Similarly, as e 4 0, e _1 [P {r c < t + e} - P {r c < t}] has the same limit. It follows that 
P {r c < •} is differentiable everywhere in the open set {t G G : c(t) > 0}, and hence its 



derivative —c'(t)gt(c(t)) + ij) c {t){t) ls t ne p.d.f. of r c on the set [34|, Theorem 7.21]. Then by 
Corollary EjZJ 

P ^ €dt > = -c>(t)g t (c(t)) + P{Tc€d *' g ^ >c ^>, 
dt dt 

which yields d2H]). □ 
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We need one more lemma before getting the second main result of this section. 

Lemma 3.6. Let S satisfy the continuous density condition. If c is continuous and non- 
increasing on (0, oo) with c(0+) > 0, then P {r c £ ^4} = for any A C (0, oo) with £(A) = 
c(A) = 0, where c(A) is the Riemann-Stieltjes integral \ J 1 {x G A n (0, oo)} dc(x)|. 

Theorem 3.7. Let c be an absolutely continuous non-increasing function on (0, oo) with 
c(0+) > 0. Suppose c is differ entiable on (0, oo) \ F for some closed set F with i(F) = 0. 
Then under the continuous density condition on S, w.p. 1, t c G (0, oo) \ F and for u G 
[0,oo) rf andv 6 [0,oc) d \{0}, 

P{X(t c -) edu, A x {r c ) Edv\r c } = Z(T e )- 1 u 1 (du,dv\T e ), (3.15) 
P {X(r c -) G du, A x (t c ) = I t c } = Z(T c y 1 v 2 (du | r c ), (3.16) 

where for t G (0, oo) \ F , 

ux(du, dv\t) = l{0< c(t) - \\u\\ < \\v\\} P{X(t) e du} U(dv), 
u 2 (du 1 1) = -c'(t)g t (c(t)) P{X(t) € du | S(t) = c(t)}, 

f c(t) _ 

Z(t) = -c'(t)g t (c(t)) + / U s (c(t) - s)P{S(t) G ds}. 

Jo 

Proof. Since c is absolutely continuous, c(F) = 0, so by Lemma I3UI r c E (0, oo) \ F w.p. 1. 
By Theorem O and Proposition E3 for t G (0, oo) \ F, u G [0, oo) d , v G [0, oo) d \ {0}, 

P {r c G dt, X{t c -) G du, A x (r c ) G dv] 

= 1 {0 < c(t) - |H| < \\v\\}dt P {X(t) G du} n(d«) = dt i^(du, dv | 
P {r c G dt, X(r c -) G du, A x (r c ) = 0} 

= -c'{t)g t (c(t)) dt P{X(t) G du | S(t) = c(t)} = dt u 2 (du 1 1). 

Integrate over u and v to get P {r c G dt} = Z(t) dt. Then f|3. 15|) follows. □ 

3.2 Results for general Levy processes with bounded variation 

For a process X taking values in M, denote 

X(t) = sup{X(s) :s<t}, X(t) = inf{X(s) : s < t}. (3.17) 
The following results will be used to validate the extension described in Section 11.21 

Proposition 3.8. Let X be a Levy process taking values in R with bounded variation and 
non-positive drift. Suppose X is not compound Poisson. Then for any a > 0, 

P {3t > such that X(s) <a all s < t, X(t-) = a or X(t) = a} = 0. 

Proof. For each t > 0, denote At = {X(s) < a all s < t}. We first show 

P {3t > s.t. X(t-) = a, A x (t) / 0} 

= = P {3t > s.t. A t , X(t) = a, A x (t) + 0} . (3.18) 
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Let LT be the Levy measure of X. Given e > 0, 

l{3t>0s.t. X(t-) = a, \A x (t)\ >e}< ^ l{X(t-) = a}. 

t:\A x (t)\>e 

Take expectation on both side and apply the compensation formula to get 
P {3t > s.t. X(t-) = a, \A x (t)\ > e} 

<II(R\ (-£,£)) / P = a} dt < n(R\ (-e,e))l7({a}), 

where [/(•) = J* P{X(t) € -}di is the potential measure of X. Since U is diffuse 0, 
Proposition 1.15] and II(R\ (—£, e)) < oo, the left hand side is for any e > 0, showing the 
first half of (|3.18p . On the other hand, since At implies X(t—) < a, 

P {3t > s.t. A t , Jf (f) = a, A x {t) ± 0} < P {3t > s.t. A^(t) = a - X(i-) > 0} , 

which is by the argument for Proposition III. 2 in Q|. This shows the second half of (|3. 18[) . 
To complete the proof, it only remains to show 

P {3t > s.t. A t and X(t-) = X(t) = a} = 0, 

or equivalently, P {r* < oo} = 0, where r* = inf{t > : At and X(t—) = X(t) = a}. 
Let T a = mf{t > : X(t) > a}. Clearly, {r* < r a } C {r* < oo} C It* < r a }. Since 
{t* = r a < oo} C {X creeps across a at r a }, has probability [cf. y, Exercise VI. 9], 
P {r* < oo} = P {r* < T a }. Let 77 ~ Exp(l) be independent of X. If P {r* < r a } > 0, then 
P {t* < t] < r a } > and hence P{X(ry) = a} > 0. However, from the fluctuation identity 
0, Theorem VI.5], Z(r?) is either constant or infinitely divisible with Levy measure 
v(dx) = 1 {x > 0} / °° i _1 e~*P {X(t) G dx} dt. In the latter case, since is diffuse, v is 
also diffuse, implying the distribution of X{rj) is continuous on (0, 00) [cf. [3^, Remark 27.3]. 
As a result, P{X(rj) = a} = 0. The contradiction implies P {r* < 00} = 0. □ 

Applying the result to X and —X respectively and using union-sum inequality, we get 

Corollary 3.9. Let X be a Levy process taking values R with bounded variation and no 
drift. Suppose X is not compound Poisson. Then for any a > 0, b > 0, 

P {3t > such that - b< X(s) < X(s) <a alls <t, X(t-) or X(t) = -b or a} = 0. 

4 Sampling of first passage event for a subordinator 

We call a function c "regular" if it satisfies the conditions in Theorem l3.71 i.e., c is absolutely 
continuous and non-increasing on (0, 00) with c(0+) > 0, and is differentiable on (0, oo)\F, 
where F is a closed set of Lebesgue measure 0. Note that if c is regular, then for any 
constant a > 0, c A a is also regular. 

Let Z be a subordinator with Levy measure II and no drift, such that 1) 11(0, 00) = 00 
and II can be decomposed as 

n(dx) = e~ qx l{x < r}A(dx) + x(dx), q > 0, r > 0, (4.1) 
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with x(0, oo) < oo, and 2) letting S be a subordinator with Levy measure A and no drift, 
its passage event across any regular function can be sampled. 

To utilize S to sample the first passage event of Z across a regular boundary c, let 
X\, X2, A3, and Q be independent subordinators with no drift, and with Levy mea- 
sures e~ qx l {x < r} A(dx), (1 — e~ qx )l {x < r} A(dx), 1 {x > r} A(dx), and %, respectively. 
Among the four, only X\ is not compound Poisson. Represent Z and S as 

Z = Xt + Q, S =\\X\\=X 1 +X 2 + X 3 withX = (X 1 ,X 2 ,X 3 )'. 

Denote rf = inf{t > : Z(t) > c(t)}, rf = inf{t > : S(t) > c(t)}, and A Q the jump 
process of Q. Table [T] describes a general procedure to sample the first passage event of Z 
across c. It actually does a little more. Given a terminal point < K < 00, it samples 
(r, Z(t— ), Az(t)), where r = A K . In particular, if c = 00 and K = 1, the procedure 
samples an infinitely divisible random variable with Levy measure II. 

Table 1: Sampling of (r, Z(t— ), A^(r)), where t = A K, c is a regular function or 00, 
< K < 00 (finite if c = 00) 

* Set T = H = D = 0, A = K, b(-) = c(-) 

1. If I? = 0, then sample (D,J) ~ (tq A j4, Aq(tq A A)), where t q = inf{t : Aq(£) > 0}. 

2. Sample ti ~ r^ Ar and set t = t\ A D. 

3. If £ = ii < D, then sample (s, v ) ~ (S(t—), As(t)) conditional on r^ r = t. 

4. If t = < t±, then sample s ~ S(t) conditional on S(t) < b(t) A r and set v = 0. 

5. Sample x ~ A"i(i) conditional on Xi(t) + X2(t) = s. 

6. If f > 0, then sample U ~ Unif(0, 1) and reset v <— vl {v < r, U < e~ qv }. 

7. Update T <^T + t. Set A = v + 1 {t = D} J, z = x + A, and update H ^ H + z. 

8. If z < b(t) and t < A, then update A A - t, D ^ D - t, &(•) «- &(• + t) - z, and 
go back to step 1 to start a new iteration; else output (T, H — A, A) and stop. 



Theorem 4.1. If c is a regular function and < K < 00, then the procedure in Table{]\ 
stops w.p. 1, and its output is a sample value of (r, Z(t— ), Az(t)), where t = rf A K. 
The claim is still true if c = 00 and K < 00. 

Proof We only prove the claim for the case where c is a regular function. The case where 
c = 00 and K < 00 can be similarly proved. 

Consider the first iteration. With A = K,D = t Q AK. Note Z(t) = X^t) for t < D. 
In step 2, with b = c, t\ is a sample value of r^. and t that of r* := r ( ^ Xr A tq A K. By 
independence, t\ ^ D w.p. 1. If t x > L>, then w.p. 1, S(D-) = S(D) < c(D) A r. Thus 
the pair (s, v) generated by steps 3 and 4 is a sample value of (5(r*— ), As(r*)) conditional 
on t* = t. Given (r*, S(t*— ), As(r*)) = (i, s, u), steps 5 and 6 sample Xl(t*— ) and 
Ai(t*) from their joint conditional distribution, where Ai is the jump process of X±. If 
t = t\ < D, i.e. S crosses c A r before D, then by part 2) of Theorem 13.11 Xi(t*—) and 
Ai(t*) are independent under the conditional distribution, following the distribution of 



13 



X\(t) conditional on S(t) = s and that of Ai(t) conditional on A$(t) = v, respectively. 
This is still true if t = D < t 1( as X(D-) = X(D) and Ai(D) = A S (D) = w.p. 1. 
By s < c(t) A r < r, P{Xi(t) € • | S(t) = s} = P{A"i(t) £ -|Xi(t) + X 2 (t) = a}, hence 
the sampling of x in step 5. Clearly, As(t) = implies Ai(t) = 0. Suppose As(t) = 
v > 0. The support of II x is within {(xi, X2> £3) : X{ > 0, at most one is nonzero}, 
such that for y > 0, n x (dy x {0} x {0}) = e _ ™l {y < r}A(dy), n x ({0} x dy x {0}) = 
(l-e-»)l{y < r}A(dy), and n x ({0} x {0} xdy) = 1 {y > r} A(dy). Then by Theorem EH 
P{Ai(t) e dy|A 5 (t) = w} =n„(dyx{0}x{0}) = 1 {y = v < r}e~ qv , hence the updating of 
v in step 6. Put together, the triplet (t, x, v) generated by the end of step 6 is a sample value 
of (r*,Xi(r*-), Ai(r*)), and A in step 7 is a sample value of A z (r*) = Ai(r*) + Aq(t*) 
and z that of Z(r*) = X±(t*—) + A^(r*). If the condition of termination is not satisfied, 
we can renew the sampling by shifting the origin to (t, Z(t)). This justifies the updating of 
A and b in step 8. Note that D is the distance in time to the current jump of Q. Once D 
becomes 0, the next jump of Q will be sampled. 

Let To = 0, and for n > 1, (T n ,H n ,A n ) the value of (T,H,A) obtained by the end 
of the n th iteration. By induction, we can make the following conclusion. For n > 1, if 
Z(T„_i) < c(T n _i) and T„_i < if, then 

T n = mi{t > T n _! : 5(t) - S , (T n _ 1 ) > [c(t) - Z(T n _i)] A r or A Q (t) > 0} A K, (4.2) 

il n = Z(T n ) and A n = A z (T n ). Evidently, Z(T n -) = H n - A n . 

To show that the procedure stops w.p. 1 and returns (r, Z{r— ), Az(t)), it suffices to 
show P{T n = t eventually} = 1. It is clear that To < r. For n > 1, if T n _i < r, then, 
since Z is strictly increasing w.p. 1, Z(T n _i) < Z{r—) < c(r) < c(T n _i). Then by (14. 2p . 
T n > T„_i. Observe that in this case, for any t G (T n _i,T n ), 

Zit) - Z(T n _i) = Xx(t) - Xi(T n _i) < 5(t) - 5(T n _i) < c(t) - Z{T n _ x ), 

giving Z(t) < c(t) and hence T ra < r. Therefore, if T n ^ r for all n > 1, then T n is 
strictly increasing and strictly less than r. Let 9 = limT n . Then 9 < r < oo. For 
n>l, the compound Poisson processes X3 and Q make no jumps in (T n ,9), and so 
5(T n+ i) - S(T n ) = X!(T n+1 ) - Xi(T n ) = Z(T n+1 ) - Z(T n ). Meanwhile, since Z(T n ) < 
c(T n ) < c(Ti) < 00, for n > 1, < Z(T n+1 ) - Z(T n ) < r. Then we get 

r > Z(T n+l ) - Z{T n ) = S(T n+l ) - S(T n ) > [c(T n+1 ) - Z{T n )\ A r. 

It is easy to see that the inequalities imply Z(T n+ i) > c(T n+ i). The contradiction shows 
that w.p. 1, T n = t for some n. □ 

5 Extensions to Levy processes with bounded variation 

5.1 Non-positive drift, positive constant level 

Let Z be a Levy process with non-positive drift, such that its Levy measure II satisfies 

poo 

n(M) = oo, / (\x\ A l)n(da;) < 00. (5.1) 
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Given constants a > and < K < oo, let = inf{i > : Z(t) > a} and r = rf A K. 
Decompose Z as Z + — Z~ , where Z^ are independent subordinators with Levy measures 

n+(dx) = i{x > o}n(dx), n~(dx) = i{x > o}n(-dx), (5.2) 

respectively, with Z + having no drift. Since the drift of Z is non-positive, if t£ < oo, then 
w.p. 1, Z makes a positive jump at t% [cf. 0, Exercise VI. 9]. Meanwhile, Z makes no jump 
at K w.p. 1. Therefore the only possible jump that Z can make at r is positive, giving 
A z (t) = A z +(t), Z+(t) = Z+(t-) + A z (t), and Z~(t-) = Z~(t). We thus consider the 
sampling of (r, Z + (t— ), Z~(t), A z (t)). Table [2] describes a procedure to do this, which 
essentially follows the description in Section 11.21 but allows a terminal point K < oo to be 
included. Note that by assumption, Z^ cannot be both compound Poisson. If II + (resp. 
II~) can be decomposed as in (|4. 1 1) . then the procedure in Table Q] can be directly called 
in step 1 (resp. 2) in Table (2) On the other hand, if one of Z^ is compound Poisson, the 
corresponding step is straightforward. 

Table 2: Sampling of (r, Z+(t-), Z~(t), A z (t)), where r = rf AIT, a is a positive constant, 
and < K < oo 



* Set T = H+ = H~ = 0, A = K, b = a 

1. Sample (t,z + ,v) ~ (r+, Z + (t + -), A z +(t + )), where r+ = r fe z+ A A. Set x = z + + v. 

2. Sample z~ ~ Z~(t). 

3. Update T <- T + t, H + <- ff + + x, ff" «- + z". 

4. If x — z~ < 6 and t < A, then update A <— ^4 — i, b ^— 6 + z~ — x, go back to step 1 
to start a new independent iteration; else output (T, H + — v, H~,v) and stop. 



Proposition 5.1. Suppose lim^oo Z(t) = oo w.p. 1 or K < oo. Then the procedure in 
Tabled stops w.p. 1, and its output is a sample value of (r, Z + (t—), Z~ (t), Az{t)). 

Proof. Let To = 0, Hq = H~ = 0, and for n > 1, let (T n , H+ , H~ ,v n ) be the value of 
(T, H + , H~ ,v) at the end of the n th iteration. By induction, for n > 1, the procedure has 
to continue into the n th iteration if and only if Z(Tfc) < a and < K for all < & < n, 
and in this case, 

T n = wi{t > T n _ x : Z + {t) - Z + (T n ^) >a- Z(T n ^)} A K > r n _i, (5.3) 
and H+ = Z + {T n ), H~ = Z~(T n ), v n = A z (T n ). By Z+(T n -) - Z+(T n _i) < a -Z(T n _i), 

Z+(T n -)-Z-(T n _ 1 )<a. (5.4) 
We show that for n > 1, if the procedure has to continue into the n th iteration, then 

sup{Z(t) : t < T n } < a w.p. 1 (5-5) 
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Consider n = 1. Since Z is right-continuous, there is e > 0, such that Z{e) < a, where Z is 
defined as in (|3.17p . Since at least one of Z^ is strictly increasing, by (|5,4p , for e < s < t < 
7i, Z(s) < Z + (t) - Z~(e) < Z + (T n -) < a. As a result, Z(t) < a. If fl53]) is not true, then 
there must be Z(T\ — ) = a. However, by Proposition 13.81 the probability for this to happen 
0. We then get (|5.5p for n = 1. For n > 2, by renewal argument, if the procedure has to 
continue into the n th iteration, then sup{Z(£) — Z(T re _i) : T n _i < t < T n } < a — Z(T n _i), 
which together with induction yields (|5.5|) . 

By assumption, r < oo w.p. 1. To finish the proof, it suffices to show w.p. 1, T n = r 
eventually The compliment of the event has two cases. The first one is that the procedure 
stops at the end of an iteration with T n ^ r. In this case, by (|5.5[) . T n < r < K. Now 
T n < t implies Z(T n ) < a, while T n < K together with the stopping rule of the procedure 
implies Z(T n ) > a. Thus Z(T n ) = a. Also by (|5.5p . Z(t) < a for all t < T n . However, 
by Proposition 13. 8\ the probability of this case is 0. The second case is that the procedure 
goes on forever. In this case, by (|5.5p . Z(T n ) < a and T n < r for all n. Then by (|5.3D . 
T n is strictly increasing with a limit 9 < r. Now for any t < 8, Z{t) < a. Meanwhile, by 
Z + (T n+ i) — Z + (T n ) > a — Z(T n ) > 0, letting n — > oo yields Z{9—) = a. By Proposition 
13.81 the probability for such 9 to exist is also 0. □ 

5.2 No drift, first exit out of an interval 

Now consider the sampling of the first exit from an interval. Let Z be a Levy process 
with no drift, such that its Levy measure n satisfies (|5.ip . Given constants > and 
< K < oo, denote I = [-a~,a + ] and let rf = inf{i > : Z(t) & I}, r = rf A K. As 
in Section [5.11 write Z = Z + — Z~ , where Z^ are independent subordinators with Levy 
measures n 1 * 1 , respectively, and no drift. Then Z makes a positive jump if it first exits / 
at o + , and a negative jump if it first exists I at —a~. We thus consider the sampling of 
(t,Z+(t-),Z-(t-),A z+ (t),A z -(t)). 

Suppose that Tl^ can be decomposed as in (|4.ip . To be specific, for a £ {±}, 

n CT (dx) = exp(-q a x)l {x < r a } A a (dx) + x CT (dx), 

where g 17 > and < r a < oo are constants and x CT (0, oo) < oo, such that, letting S a be 
a subordinator with Levy measure A CT and no drift, its passage event across any positive 
constant level can be sampled. In the case where Z a is compound Poisson, we simply set 
S a = and the time of the first passage of S a across any positive boundary to be oo. Note 
that, since n(M) = oo, at most one Z a is compound Poisson. 

To utilize to sample (r, Z + (r— ), Z~(t— ), A z + (r), A z - (r)), the idea is similar to 
that in Section UJ Table [3] describes a procedure to do this. In each iteration, we have to 
monitor two passage times, i.e., the times when 5°" cross W A r CT , a € {±}, respectively, 
where b a is a constant obtained from a° ' . To simplify notation, denote by T^ Ar (a) the first 
passage time of S a across b a A r a . Represent S a = \\X a \\ = X? + X% + X$, where X a = 
{X° , X% , A~3 )' and Xf are subordinators with Levy measures exp(— q a x)l {x < r a } A°"(d2;), 
[1 - exp(-g ff x)]l {x < r a } A' J (dx), and 1 {x > r a }A a (dx), respectively. All of X? , i < 3, 
a € {±} are assumed to be independent with no drift. 
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Table 3: Sampling of (r, Z+(r-), Z~(t-), A z + (t), A z - (t)), where r = rf A K, I = 
[—a~, a + ], a~, a + > are constants, and < K < oo 



* Set T = H + = R~ = D = 0, A = K, b + = a + , b~ = a". Let Q be a compound 
Poisson process with Levy measure x + + X~ ■ 

1. If £> = 0, then sample (D, J) ~ (tqAA, A q (t q AA)) and set J+ = JVO, J" = (-J)VO, 
where r Q = wf{t : A Q (t) ^ 0}. 

2. For <T G {±}, sample t CT ~ T^ Ar (a). Set t = t + At~ AD. (Note: w.p. 1, £± and I? arc 
different from each other.) 

3. For a G {±}, sample (x a ,v CT ) ~ (Xf (£-), A X f (*)) conditional on T+ArM) = t, 
by applying steps 3-6 in Table [T] to X a . 

4. Update T <- T + t. For a G {±}, set A a = v a + 1 {t = D} J CT , z a = x a + A CT , and 
update H a <- H a + z u . 

5. If z + - z~ G (-b~,b + ) and t < A, then update A 4- A - t, D D - t, b + ^ 
b + + z~ — z + , b~ <— b~ + z + — z~ , and go back to step 1 to start a new iteration; else 
output (T, - A + , - A", A+, A") and stop. 



Proposition 5.2. Suppose Z =£ 0. Then the procedure in Table stops w.p. 1, and its 

output is a sample value of (r, Z + {t— ), Z~{t— ), A z +(t), A^-(r)). 

Proof. First, r < oo w.p. 1 as lim^oo |Z(i)| = oo [1, Theorem VI. 12]. Let To = 0, Hq = 
Hq = 0, and for n > 1, let (T n , H+, H~, A+, A~) be the value of (T, H + , H~ , A + , A~) 
obtained by the end of the n th iteration. By induction and the same argument as in the 
proof of Theorem 14. 1| for n > 1 , the procedure has to continue into the n th iteration if and 
only if Z(Tk) G (—a - , a + ) and < if for < k < n, and in this case, 

T n = inf{t > T n _! : 5+(t) - 5+ > a+ - Z(T n _i), 

or S~(t) - S^(T„_i) > a~ + Z(T n _i), or A Q (i) > 0} A if > T n _i, 
H+ = Z + {T n ), H- = 2T(T n ), A+ = A z+ (T n ), A~ = A z -(T n ), 

and, as in the proof of Proposition 15. 1( Z + (T n — ) — Z~(T n _i) < a + , sup{Z(t) : t < T n } < 
a + , and likewise, by considering —Z(t), Z~(T n —) — Z + {T n ^\) < a - , inf{Z(t) : t < T n } > 
—a~. The rest of the proof follows that of Proposition l5.lt except Corollary 13.91 is used. □ 

6 Sampling issues involved 

The procedures in previous sections involve several types of sampling, some of which are 
standard, while the others have to be dealt with on a case-by-case basis. Consider the 
procedure in Table [TJ The main task of its step 1 is 

sample the first jump of a compound Poisson process on (0, oo). (6-1) 
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The rest of the procedure requires a subordinator S with infinite Levy measure A and ana- 
lytically tractable properties be available. Under this prerequisite, given regular boundary 
c, t > 0, and < s < r, the main tasks of steps 2-5 are to sample 

rf = inf{t > : S(t) > c(t)}, (6.2) 

(S(t—), As(t)), conditional on = t, (6-3) 

S(t), conditional on S(t) < c(t), and (6-4) 

conditional on X^t) + X 2 (t) = s, (6.5) 

respectively, where X±, X2 are independent subordinators with no drift and with Levy 
measures 1 {x < r} e~ qx A(dx) and 1 {x < r} (1 — e~ qx )A(dx), respectively. The procedures 
in Tables [2] and [3] also boil down to (I6.ip - (l6.5p . 

For (|6.ip . recall that if a compound Poisson process Q has Levy measure x^O, then the 
time and size of its first jump are independent with p.d.f. qe~ qt l {t > 0} and distribution 
x/q, respectively, where q = J dx [cf. H, [ifl. For complicated \i a rejection sampling 



method known as "thinning" can be used [cf. [9|, [l6|. Let /ibea Levy measure such that 1) 
a = f dfi < 00 is readily available, 2) a fi is easy to sample, and 3) dx = gdfi for some 
function g < 1 that is easy to compute. Then the thinning based on /1 is as follows. 



* 


Set t = 




1. 


Sample s from the distribution with p.d.f. ae~ as l {s > 0}. Update t 


<- t + s. 


2. 


Sample x from the distribution a" 1 ^ and U ~ Unif(0, 1). 




3. 


If U < g(x), then stop and output (t,x) as a sample of the time and 


size of the first 




jump of Q; else go back to step 1. 





For (|6,2p . since S is strictly increasing w.p. 1 and c is non-increasing, 

P{r c 5 <t} = P{S(t) > c(t)}, (6.6) 



with each side being continuous and strictly increasing in t > 0. If the probability distri- 
bution of S(t) is analytically available, then may be sampled by inversion method, i.e., 
sample U ~ Unif(0, 1) and return the unique value of t satisfying P{S(t) > c(t)} = U [cf. 0]. 
Alternatively, if 5 has scaling property, it can be utilized to sample r^f. Both possibilities 
are illustrated later. The sampling for (|6.3p heavily relies on the results in Section [3j and its 
detail need be dealt with on a case-by-case basis. The sampling for (|6.4p has the following 
generic solution: keep sampling x ~ S(t) until x < a. However, by utilizing the structure 
of S(t), it is possible to make the sampling significantly more efficient. 

Finally, some comment on (|6.5p . Give t > 0, for the non-trivial case q > 0, if S(t) has a 
bounded p.d.f. gt, then in principle rejection sampling can be used 0, 0, [H, H2] . Indeed, 
since the Levy measure of X\(t) is e~ qx v{dx), where v(dx) = tl{x < r}A(dx) is that of 
Xi{t) + X 2 (t), P{Xi{t) G dx} oc e~ qx P{Xi(t) + X 2 {t) G dx}. Recall 5 = X x + AT 2 + X 3 , 
where X3 has Levy measure 1 {x > r} A(dx). Since X^(t) is either or > r, X\{t) + X 2 {t) 
has p.d.f. g t (x)/P{X 3 (t) = 0} on (0,r]. It follows that Xi(t) has a p.d.f. on (0,r] in 
proportion to e~ qx gt(x), and given s G (0,r], 

P{Xi(t) edx\X 1 (t) + X 2 {t) =s}(xe~ qx g t (x)P{X 2 (t) Ss-dx}. (6.7) 
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Thus to sample (|6.5p . we may keep sampling v ~ X^ify an d U ~ Unif(0, 1) until v < s and 
sup x gt(x) ■ U < e~ q ( s ~ v "> gt(s — v) and then output s — v. Here, since ^(t) is compound 
Poisson, its sampling is standard 0, [3]. However, a problem is that ^ can be hard to 
evaluate. To get around the problem, the structure of S(t) needs to be exploited. 



7 Exponentially tilted upper truncated stable Levy density 

Consider the measure n(dx) = 1 {x < r} e~ gx A(dx) + x(dx) specified in (jl.ip . where 
A(dx) = 1 {x > 0}7x _1_a dx with a € (0,1). Let c be a regular function as defined in 
Section [H As an application of the procedure in Table [H we next show an algorithm 
to sample the first passage event of a subordinator Z with Levy measure n and no drift 
across c. Using the procedure in Table [2] or [31 an algorithm can be similarly devised to 
sample the first passage event of a Levy process with non-positive or no drift across a 
constant level or interval, when its Levy measure is 1 {0 < x < r + }^/ + e~ q x x~ 1 ~ a dx + 
1 {— < x < 0}^f~e~ q x \x\~ 1 ~ a dx + x(dx) with a^ 1 € (0, 1). We omit detail on this. 

Let S be a stable subordinator with Levy measure A, so that for A, t > 0, E[e~ xs ^} = 
exp{— t'jTil — a)a~ 1 X° 1 }. By scaling of time, assume 7 = a/r(l — a) without loss of 
generality. Then S(l) is a "standard" stable variable with p.d.f. 

fix) = t r— / h(x, 6) d6, where 

(1 - a)TT Jo 

h(x, 6) = 1 {x > 0} hoiff) x- 1 ^ 1 -^ exp{-h (9) x~ a/{1 - a) }, with ( 7A ) 

h (0) = sin[(l - a)^][sin(a0)] a /( 1 - Q )(sin0)- 1 /( 1 - Q ). 

The sampling of S(l) is well-known 0,0,113]. Define function 

if>(x) = l{x + 0}x _1 (l - e~ x ) + l{x = 0} . 

Let < K < 00 and r = inf{i > : Z(t) > c(t)} A K. An algorithm to sample 
(r, Z(t— ), Az(t)) is as follows. 

* Set T = H = D = 0, A = K, &(•) = c(-), M a = (1 - a) l ~ l l a a~ l ~ l / a e~ l / a . 

1. Sample (D,J) as in step 1 in Tabled! 

2. Sample 5(1). Set t x such that t\ /a S(l) = b(h) A r, t = h A D, and 2 = &(t) A r. 

3. If i = h < D, then set 



^0 



d(6(«)Ar) -fz 1 -* 



du 



u=t «(! - ") 



and do the following steps. (Note: w.p. 1, b{u) A r is differentiable at t with a non- 
positive derivative.) 

(a) Sample 1 € {0, 1} such that P{l = i} = Wi/{wq + w\). If t = 0, then set s = z, 
v = 0; else sample /3i ~ Beta(l,l — a), fa ~ Beta(a, 1), and set s = /3iz, 
u = (z - s)//3 2 . 
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(b) Sample ■& ~ Unif(0,7r) and U ~ Unif(0,l). If U > h((t)~ 1 / a s,^)/M a , then go 
back to step 3(a). 

4. If t = D < ti, then sample <S(i) conditional on S(i) < z. Set s = 5(t) and t> = 0. 

5. Set 



k\T(l + fc(l - a)) 



k=0 



Then do the following steps. 

(a) Sample k G {0, 1, 2 . . .}, such that P{k = k} = C^/C. 

(b) If k = 0, set a; = s, q = 1; else sample /3 ~ Beta(l, k(l — a)) and (cji, . . . , oj k ) ~ 
Di(l — a, . . . , 1 — a), and set x = s/3, g = Y\ i ip(q(s — x)u>i). 

(c) Sample ~ Unif(0, vr) and 17 ~ Unif(0, 1). If M a U > ge~ qx h{t~^ a x, 0), then go 
back to step 4(a); else output x to the next step as a random sample of X\(t—) 
conditional on X\(t—) + X2(t—) = s. 

6. The rest is the same as steps 6-8 in Table [TJ 

We next justify the algorithm. All its steps correspond 1-to-l to those in Table [H 
and only steps 2, 3 and 5 contain new detail. Denote a = b A r, which is regular. With 
t\ being the unique solution to t 1/ ' Q S'(l) = a(t), from (I6.6D and scaling property of 5, 
P{ T( f < t} = P{t^ a S(l) > a(t)} = P{h < t}, leading to step 2. Given rf = t x and 
r* = t, where r* = A D, by step 3 in Table [H if t = t\ < D, then we need to sample 
(S(t— ), As(t)) conditional on = t. From Theorem 13.71 

P{S(t-) G ds,A s (t) G dv | rf = t} oc # t (s) [w <5(ds - z)5(dv) + w±p(s,v) dsdv], 

with z = a(t), t^o = |«'(i)| 5 w\ = 7z 1_a /[a(l — a)], and p is the following p.d.f. 

p(s, v) = 1 {0 < z - s < v} a{\ - a)z~ 1+a v~ l - a . (7.2) 

Define random vector (i,(,V), such that P{l = 0} = 1 — P{l = 1} = wo/{wq + wi), 
P{( = 2 , V = 1 1 = 0} = 1, P{C € ds, V G dv | l = 1} = p(s, i>) ds dv. Let ~ Unif(0, vr) 
be independent of (l, £, V). Then by <?t(s) = i _1//o /(i _1 / Q s), where / given in (17. 



P{S(t-) G ds, A s (t) G dv | rf = t} 

h(t~ 1/a s, 9)P{tf G d9, l G di, C G ds, V G dv}, 



ex 



(7.3) 



with the integral only over 9 and i. It is seen that £ ~ (1 — JTj 1 ^ 1 Q ^)z and conditional on 

C, V ~ (z- C)U2 l/a , with C/i, U 2 i.i.d. ~ Unif(0, 1). Thus step 3(a) samples (0,t,£V). 
Next, for x > and G (0, 7r), by change of variable s = x _Q 7( 1_Q ) in the expression of h, 



h(x, 9) < sup 

6>e(0,7r) 



h {9) x sup(s l/a e- ho{e)s ) 

s>0 



a 



-l/a.-l/a 



sup /lo(#) 

6»G(0,7r) 



1-1/q 
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Since sin(i0)/sin(6>) > t for 9 G (0, vr) and t G (0,1), h (6) > (1 - aOa ^ 1 " ), giving 
h(x,0) < M a . Thus, step 3 is a rejection sampling procedure for the distribution which 
is proportional to h^^s, 6»)P{i9 G d(9, z G di,C G ds, V G dv}. Then by £L3|), is a 

sample of {S(t— ), Ag(t)) conditional on r„ = t. The entire step 3 is now justified. 

By step 5 in Table[[J given (r*, 5(r*— )) = (i, s) with s G (0, r], we need to sample X\(t) 
conditional on X\(t) + X2(t) = s, where X\ and X2 are independent subordinators with 
Levy measures 1 {x < r} e~ qx A(dx) and 1 {x < r} (1 — e~ qx )A(dx), respectively. By (|6.7|) . 

P{X x (t) G dx I X^t) + X 2 (t) = s} oc e- 9a; ^(x)P{X2(t) G s - dx}, < s < s. 

Since X2(t) is compound Poisson with Levy density A(x) = 1 {0 < x < r} t^(\—e~ qx )x~~ 1 ~ a , 



P {X 2 (t) G s - dx} oc 1 {0 < x < s} 



<5(s - dx) + ^ 



fe=i 



A* fc (s - x) dx 
fc! 



where X* k is the /c-fold convolution of A. For w > and A; > 1, 



\* k (w) = w k ~ l [ Y[\(wv t )a k (d 

^ i=i 



where is the measure specified in Section [2l Since < w < s < r, by the definition of ip 
and Dirichlet distribution, for any v = (v%, . . . ,v k ) with i>j > and = 1, 



nA(^)=(n) fc n 



i=i 



1 _ e -qwvi 
\ (wvi) 1+a 



1 



(t 1 ) k q k w- ka J{^qwv l )W- 



i=l 



i=l 
k 



W 



-ka 



T(k(l - a)) ax 



where denotes the p.d.f. of Di(ai, . . . , a^), with all Oj = 1 — a. Denote by u k 
(oJki, ■ ■ ■ ,u>kk) a random variable following the Dirichlet distribution. Then 

[t iq r(i-a)] k l 



A* fc (s - x) dx = {s - x)^ 1 "")- 1 dx 



r(fc(l-a)) 



JV(?( 



S - XjWfejj 



i=l 



Note that 1 {0 < x < s} fc(l - a){s - x )H^-»)-i / s k(i-a) ig the p d f of s ^ where ^ ^ 
Beta(l, k(l — a)). Then, with C k the same as in the algorithm, we get 

k 



X* k {s - x) dx = k\C k P{sp k G dx}E 



~[tp(q(s - x)uj ki ) 



The above identity also holds for k = 1. Combining with (|7,ip . we get 
P{Xi(t) Gdx|Xi(t)+X 2 (t) = s} 



oce qx g t (x) 



5(s - dx) + ]T C fe P{s/3fc G dx}E 



fc=i 



J^(<?(s - x) 



Li=l 



oc / e-^7i(t- 1/a x,0)d0i5(s-dx) + ^C fc P{s/3 fc Gdx}E 



fc=i 



n^( 



s - x)u ki ) 
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Starting from the last integral expansion, the treatment is similar to step 3. Define random 
vector (re, £, lj), such that re G {0,1,2,...} with P{re = k} oc C^, conditional on k = 0, 
C = s, oj = 0, and conditional on re = /e > 1, £ ~ s/3fc and to ^ u)k are independent. For 
any x G [0, s] and i« G U^iIS^, define g(x,w) = Ili=i ^(<z( s ~ x ) w i)i with k equal to the 
dimension of w. Finally, let t? ~ Unif(0,7r) be independent from (k,(,ui). Then 

P{Xx(t) edx\X 1 (t) + X 2 (t) = s} 

oc J e~ qx h(t~ l / a x, 9)g(x, w)P{i) G d9, re G dk, ( G dx, u) G du;}, 

where the integral is only over 0, k, and Based on this, it is seen step 5 is a rejection 
sampling procedure of Xi(t) conditional on Xi(t) + X^it) = s. 

8 Finite mixture of exponentially tilted upper truncated sta- 
ble Levy densities 

Given / > 2, r, 71, . . . , 71 > 0, q > 0, and ai, . . . , a\ G (0, 1), let 

U(dx) = ^2if i{x)dx + x(dx) with ipi(x) = 1 {0 < x < r} j ie ~ qx x' 1-0 ' . (8.1) 

The seemingly more general case where different (fi(x) have different and is actually 
covered by (|8.ip . once we set f = minrj, q = max^j, ^(x) = 1 {1 < x < f} ^fie~ qx x~ 1 ~ ai 
and x(dx) = x{^ x ) + J2i[ { Pi( x ) ~ <Pi( x )]dx. As in last section, we shall focus on the 
sampling for subordinators. The extension to Levy measures ^2 i=1 ifi(±x) dx + x(dx) is 
rather straightforward. 

To apply the procedure in Table [TJ let Xij, i = j = 1,2,3, and Q be in- 

dependent subordinators with no drift, such that Xn, Xi2, and X^ have Levy densi- 
ties tpi(x), 1 {0 < x < r}7j(l — e~ qx )x~ l ~ a \ 1 {x > r}^iX~ l ~ a \ respectively, and Q has 
Levy measure x- Let Si = Xn + JQ2 + AQ3 and Z = Yli-^n + Q- Then Si, . . . ,5/ 
are independent stable processes with Levy densities 7iX _1_a % respectively, and no drift, 
and Z a subordinator with Levy measure II and no drift. Let 5 = (5i,...,5j). Un- 
like previous sections, 5 is multidimensional. Let S = \\S\\ = S\ + • • • 4- Si. Then 
s (0 ~ Ei=i* 1/ai5 'i( 1 )- Let d i = [«i/7*r(l - a*)] 1/Qi . Then each d^l) has Laplace 
transform E[e~ Xdi ^W] = exp(— A Qi ), A > 0, whose p.d.f. fi is given in (j7.1jl . Let be hi the 
function in the integral representation of fi in (|7.ip . 

Let < K < 00 and r = inf{t > : Z(t) > c(t)} A K. An algorithm to sample 
(t, Z(t— ), Az(t)) is as follows. 

* Set T = H = D = 0, A = K, &(•) = c(-), Mi = (1 - c^-V^a- 1-1 / ^- 1 /**, i < / 

1. Sample (D, J) as in step 1 in Table [TJ 

2. Sample 5(1). Set h such that J2i=i t{ /ai Si(l) = b(h)Ar, t = t±AD, and z = b(t)Ar. 

3. If t = h < D, then set 



w 



1 d(b(u) A r) 



r(i) du 



7i^- ai T(l - ai) 



u=t ajly + 1 - ai) 
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and do the following steps. 

(a) Sample i 6 {0,1,..., /}, such that P{t = j} = Wjj (wq + w\ + ■ ■ ■ wi). Sample 
to = (lj±, . . . ,oji) ~ Di(l, . . . , 1). If i = 0, then set Si = zuii, i < I, and v = 0. 
If i > 1, then sample /3 ~ Beta(I, 1 — a t ), (3' ~ Beta(a t , 1), and set Sj = z/3u>i, 
i < I, and v = (z — si — ■ ■ ■ — Si) / /3' . 

(b) Sample ~ Unif(0, vr) and J7 ~ Unif(0, 1) independently. If U > 
llLiihiidit^/^Sutf^/Mi}, then go back to step 3(a). 

4. Ut = D <h, then sample 5i(l), . . ., 5/(1) conditional on Y,l=i £> 1/Qi £i(l) < z - Set 
s = Y.l=i -D 1/Ql Si(l) and t; = 0. 

5. For % < I, sample xn from the distribution of Xn{t—) conditional on Xn(t— ) + 
Xi2(t—) = Si, by executing step 5 in the algorithm in Section [7] with the corresponding 
parameters. Set x = x\\ + xi\ + • • • + xj\. 

6. The rest is the same as steps 6-8 in Table [TJ 

To justify the algorithm, denote a(t) = b(t) A r. From (|6.6[) and scaling property of Si, 
P{r a E <t} = P{E(t) > a(t)} = P{ELi * 1/ai -Si(l) > a(t)}, which leads to step 2. 

Given = t\ and r* = t, where r* = r E A D, denote 2 = a(t). Denote by g t the p.d.f. 
of By Theorem 13.71 for s = (s\, . . . ,sj) with Sj > 0, and u,v > 0, 

P{S(t-) € ds, E(t-) € du, A s (i) G dv | = t} 
oc g t {u) P{S{t) G ds | £(t) = u} 

|a'(t)| <S(dit - z)<y(dt;) + 1{0<z-u<v}^ 7 * 

2=1 

Since each <Si(t) has p.d.f. proportional to fi{dit~ x / ai x), 

g t {u) P{S(t) G ds | E(t) = u} oc Jl fiickt-^Si) X P{u W G ds} x 

with oj = (wi, . . . ,ujj) ~ Di(l, . . . , 1). Let u>o, u>i, ... be as in the algorithm. Then 
x \a'(t)\5(du - z) 5(dv) = z I ~ l w 5(du - z) 5{dv), 



T(I) 



1 7j du dv 



. xl{0<z-u<v} '\ = z^WiPizpi G du, {z - u)/& G dv}, i > 1, 

where ~ Beta(1, 1 — c^) and fi[ ~ Beta(aj, 1) are independent. The above identities 
together with (jT.lt) yield 

P{S{t-) G ds, E(t-) G du, A s (t) G dv | r a s = t} 

oc / TT h i (dit~ 1/ai Si, Oi) d6>i • • • d0j P{uw G ds} 

J0i£[0,Tr] fj[ 

I 

w 5(du - z) 5(dv) +^WjP{zPi G du, (z - G dv} 



i=l 
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where the integral is only over 0j's. Starting here the treatment is similar to Section [71 It 
is then seen that step 3 samples (s,v) ~ (S(t—), As(i)) conditional on = t < D. 

From step 4 in Table [T] and the scaling property of Si , it is easy to see step 4 sam- 
ples E(-D) conditional on D < t~. Finally, note the ultimate goal of step 5 in Table [T] 
is to sample ) conditional on r* = t. Since we have now sampled (s,v) ~ 

(S(t— ), Ae(£)) conditional on t* = i, it suffices to sample {Xn(t— ),..., Xn(i— )) condi- 
tional on (S(t— ), As(t), r*) = (s,v,t). Since in this case, A^i— ) + X i2 (t— ) = Sj and 
Xi3(t—) = 0, the conditional sampling is equivalent to that of X(t— ), where X = (Xij,i = 
l,...,I,j = 1,2,3). By Theorem EH if t = h < D, then 

P{X(t-) G dx I r* = t, E(t-) = ||s||, A E (t) = v} = P{X(t) G dx I E(t) = ||s||}, 
P{S(t-) G ds I r* = t, E(i-) = ||s||, A E (i) = w} = P{5(t) G ds | E(t) = ||s||}, 

yielding P{X(t-) G dx j r* = t, 5(t-) = s, A s (t) = v} = P{X(t) G dx\S(t) = s}. The 
identity still holds if t = D < t±. By independence, the right hand side is 

1 

Y[P{X a (t) G dxa,X i2 (t) G dx i2 ,X i3 (t) G dx l3 \Si(t) = sj, 

i=l 

so to sample (Xu(t-),...,X;i(f-)) conditionally, it suffices to sample -Xii(t) indepen- 
dently, each conditional on Si(t) = Sj. Therefore, step 5 in the algorithm in Section [7] can 
be used. Once Xn(t) are sampled, their sum is a sample value of ^2iXn(t). 

9 Upper truncated Gamma Levy density 

Let IT(dx) = if(x) dx + x(dx), where 

<p(x) = 1 {0 < x < r} e~ x x~ l , r > 0, (9.1) 

and x is a finite measure on (0, 00). For a Levy measure of the more general form n(dx) = 
(p{x)dx + x(dx), where <p(x) = 1 {0 < x < r} r ye~ qx x~ 1 , 7, q > 0, the sampling of the 
first passage event can be reduced to that for IT. Indeed, if Z has Levy measure II, then 
Z(t) = qZ(t/j) has Levy measure (|9.ip . with r = qr, x(dx) = x(dx/q) therein. As a result, 
the first passage event of Z across c can be deduced from that of Z across qc(t/^f). 

The sampling of the first passage event for II is somewhat simpler than those in previous 
sections, as exponential tilting is "built in" the Gamma process. Let X\, X 2 , and Q be inde- 
pendent subordinators with no drift and with Levy measures (p(x) dx, 1 {x > r} e~ x x~ 1 dx, 
and Q, respectively. Then S = X\ + X 2 is a Gamma process and Z = X\ + Q has Levy 
measure II and no drift. Let < K < 00 and r = inf{i > : Z(t) > c(t)} A K. An 
algorithm to sample (r, Z(t— ), A^(r)) is as follows. 

* Set T = H = D = 0, A = K, &(•) = c(-) 

1. Sample (D, J) as in step 1 in Table [TJ 

2. Sample U ~ Unif(0, 1). Set h such that / 6 ~ l)Ar x^e'* dx/Tfa) = U, t = t x A D, 
and z = b(t) A r. 
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3. If t = ii < D, then set 



d{b(u) A r) 

Wq = — 



du 

h\ (x, v) = 1 {0 < z — x < v < z} 



2B(t,l/2)z 1 

wi = , w 2 = - 

e t 



u=t 

^(l-x/z) 1 / 2 ln[(l -x/z)- 1 } 



2/e 

ZG 

h 2 (x, v) = 1 {0 < z — x < z < v} . 

V 

and do the following steps. 

(a) Sample i G {0, 1}, such that P{t = i} = Wi/(w$ + W\ + w 2 ). If t = 0, then set 
x = z, v = 0, 77 = 1; if i = 1, then sample j3 ~ Beta(t, 1/2), £ ~ Unif(0, 1), and 
set x = z/3, v = z(l — /3)£, 77 = h\(x,v); if 1 = 2, then sample /3 ~ Beta(i, 1), 
£ ~ Exp(l), and set x = z/3, v = z + £, r? = /i2(x, t>). 

(b) Sample C/ ~ Unif(0, 1). If £/ > r], then go back to step 3(a). 

4. If t = D < t±, then sample 7 ~ Gamma(Z), 1) conditional on 7 < z. Set x = 7, u = 0. 

5. The rest is the same as steps 6-8 in Table HJ 

To justify the algorithm, denote a(t) = b(t) A r. From (|6.6p and S(i) ~ Gamma(t, 1), 
PK 5 < t} = P{5(t) > a(t)} = ^- / x*-^- dx, 

1 W Ja(t) 

which is a continuous and strictly increasing function. The sampling in step 2 is then the 
standard inversion [cf. @|. Given = t\ and t* = i, where t* = rf A D, by step 3 in 
Tabled! if t = t\ < D, then we need to sample (S(t— ), As(/j)) conditional on 7™ = t. Let 
<7t denote the p.d.f. of Gamma(t, 1). Let z = a(t). From Theorem 13.71 for x,v > 0, 

P{S(t-)edx,A s (t)edv\r^ = t} 

e~ v 

oc \a'(t)\gt(z)5(dx — z)5(dv) + 1{0<z — x<v} gt{x) dx dv 

v 

oc \a'(t)\5(dx — z)5(dv) + q±(x,v) dxdv + q2(x,v) dxdf, 

where, letting q(x,v) = gt{x)e~ v /[vgt(z)], qi(x,v) = 1 {0 < z — x < v < z] q(x, v) and 
q 2 (x, v) = 1 {0 < z — x < z < v} q(x, v). Now q(x, v) = (x/z) t ~ 1 e z ~ x ~ v /v. Let 

/ s r , (x/zY-Ul -X/Z)" 1 / 2 1 

pi(x,v) = 1{0 < z-x < v < z}^-' ' v 



p 2 (x, v) = 1{0 < z — x < z < v} 



B(t,l/2)z v\a[{l- x/z)~ 

t(xlzf- x e z - v 



z 

For i = 1,2, pi(x,v) is a p.d.f. and qi{x,v) = Wihi(x,v)pi(x,v), where Wi and hi are defined 
in the algorithm. It is easy to check hi(x,v) < 1 for i = 1,2. Define ho(x,v) = 1. Define 
(i, C, V) such that 1 € {0, 1, 2} with P{l = i} = Wi/(wo + w\ + ^2), conditional on 1 = 0, 
C = -2 and V = 0, and conditional on z. = i € {1, 2}, (£, y) has p.d.f. pi. Then 

PIS 1 ^-) G dx, A s (t) £ dv\T% = t} = J hi(x,v)P{i £ di,( G dx, 1/ G dv}, 
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where the integral is only over i. It is easy to check that when (C, V) has p.d.f. pi, then 
(C, V) ~ (zjSt, z (! - A) 17 ), wh ere /3 t ~ Beta(i, 1/2) and 17 ~ Unif(0, 1), and when (C, V) has 
p.d.f. p2, then (£, V) ~ (z/3' t ,z + £)> where ~ Beta(i, 1) and £ ~ Exp(l). Then step 3 
in the algorithm is a rejection sampling procedure of (S(t-), As(t) conditional on r„ = t. 
Since Sft—) sampled by step 3 or 4 is exactly Xi(t—), there is no need for a counterpart of 
step 5 in Table [TJ We can directly go to steps 6-8 in Table [TJ 



Appendix 

To prove Lemmas 13.41 and 13. 6[ we start with two more lemmas. 
Lemma A.l. For t > and < a < b < oo, let 

ra rb 

L 1 ft,a,b) = [U s (a - u) - U s fb - u)]g t fu)du, L 2 (t,a,b) = U s (b - u)g t {u) du. 

JO J a 

Then for any E = [to,t±] C (0, oo) and I = [a, (3] C (0, oo), 

limsup {Li(t, a,b) : t E E, a,b £ I, < b — a < r} = 0, i = 1, 2. (A.l) 

Proof. Since II5 is non-increasing, for < b — a < r, L\(t,a,b) < Li(t,a,a + r). Fix 
e G (0, a/2) and let h(u) = [lis (a — u) — Ug(a + r — u)]g t (u). Then Lift, a,a + r) = Ji + J 2 , 
where Ji = Jift, a, r) = h, J2 = J2(t, a, r) = J £ a h. We have 

Ji< U s fa - u)g t (u) du < U s {a - e) / g t fu) du 
Jo Jo 

= n 5 (a - e)P {Sft) <e}< U s fa - e)P{Sft ) < e} 
and letting M = sup{gt(u) : t G E, e < u < f3}, 

/a _ _ rfi _ 

[U s fa -u)- U s fa - u + r)] du < M J [U s fu) - U s fu + r)] du, 

By the assumption on gtfx), M < 00. Also, Usfu) du = J °°fv A f3)Usfdv) < 00. Then 
by monotone convergence, as r J. 0, J2 — > uniformly for ft, a) G E x I. As a result, for 
Li, the limit in ([A~T]) is no greater than U s fa - e)P{Sft ) < e}. Since P{5 , (t ) > 0} = 1 
and e is arbitrary, the limit is equal to 0. Thus (jA.ip holds for L\. 

Next, fixing e G (0, a), by change of variable, for t G a, b G / with < 6 — a < r, 

L 2 ft,a,b)= I Usfdx) / g t < / n s (dx) / g t + U s fdx) / fft 

JaV(b-x) Jo Jb-x Je Ja 



JaV(b~ 

< M' 



[ xU s fdx) + rU s fe) 
Jo 



where M' = sup{gtfu) : t £ E,a — e < u < /?}. Therefore, as r \. 0, the limit for L2 in 
(jA.lj) is no greater than M' f Q £ xUsfdx). Since e is arbitrary, the limit is 0. □ 
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Lemma A. 2. Let h be a bounded function on (0, oo) x (0, oo). For a, t G (0, oo), define 
H(a, t) = j 1 {u < a < x} h(u, x) P{S(t) G du} U s (dx - u). 

Then under the continuous density condition, H is continuous on (0, oo) x (0, oo). 

Proof. Without loss of generality, suppose \h(u,x)\ < 1. It suffices to show H is continuous 
on any R = [a, 0\ x [to, h] C (0, oo) x (0, oo). Let (a, s), (b, t) G R. Then 

\H{b,t) - H(a,s)\ < \H(b,t) - H(a,t)\ + \H(a,t) - H(a, s)\ 

Let L\ and L 2 be as in Lemma lA.li Let a' = a A b and b' = a V b. Then 

\H(b,t)-H(a,t)\ < J \l{u < b < x} - l{u < a < x} | P{S(t) edu}U s (dx-u) 

< J (l {u < a' < x < b'} + 1 {a' < u < b' < x}) P {S(t) G du} U s (dx - u). 

The right hand side is Li{t,a',b') + L^it, a', b'). Then by Lemma [A.li as (b,t) — > (a,s), 
Li(t, a', b') + L 2 (t, a', b') 0, giving H(b, t) - H(a, t) 0. 

It only remains to show H(a,t) — H(a,s) — s> 0. Given e G (0,a), let M = sup{gt(u) : 
u G [a- e,/3],t G [*o,*i]}- Then 

\H(a,t) — H(a, s)\ < J 1 {u < a < x} \gt(u) — g s (u)\Tls(dx — u) du 

= J 1 {u < a} \g t (u) - g s (u)\U s (a - u) du. 

Bounding the integral on [a — e, a] and [0, a — e] separately, we obtain 

\H(a,t)-H{a,s)\<2M J U s (u) du + TT s (e) J \g t {u) - g s {u)\ du, 

Let t —¥ s. Since point-wise convergence of gt to g s implies convergence in total variation, 
]\m\H(a,t)-H(a,s)\ <2M TT s (n) du < 00. Letting e ->• gets H(a, t)-H(a, s) 0. □ 

Proof of Lemma \3.4\ Apply (13. 12H and Lemma lA.21 with h(a,t) = 1 therein. □ 

Proof of Lemma \3.6l Let G c = {t > : c(t) > 0}. Then G c is a non-empty open interval 
and P {r c G G c } = 1. To prove the lemma, it suffices to show that for any < to < t\ < 00 
with [to,tx] C G c and A C (to,h) with 1(A) = c(A) = 0, P {r c G A} = 0. Let a = c(h) 
and /3 = c(to). Given e > 0, A can be covered by at most countably many disjoint intervals 
(di,bi) C (to,ti) such that XX^« — fli) < e an d S[ c (°«) — c (^«)] < e - For each i, 
P{r c G (oi,6i)} < P{S{a t ) < c( ai ), SQh) > c(6 4 )} 

< P {cQh) < S(oi) < c( ai )} + P {5(oi) < c(6i) < 5(6,)} 
= P {c(6i) < S( ai ) < c(ai)} + P {r c(b!) G (oi, 6,)} . 

Since aj,6j G [toj*i] an d c(aj),c(6j) G [a, /3], by the continuous density condition and 
Lemma [331 the right hand side is no greater than M\[c(ai) — c(6j)] + M<i(bi — etj), where 
Mi = supgt(x) and M 2 = sup^ x (i) over (t,x) G [io,*i] x [a, Therefore, P {r c G ^4} < 
P {t c G (a,, 6j)} < (Mi + M2)e. Since £ is arbitrary, this yields the proof. □ 
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